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I Abstract 

^^ An error control technique aimed to assess the quality of smoothed 

vN finite element approximations is presented in this paper. Finite element 

?— ( techniques based on strain smoothing appeared in 2007 were shown to 

provide significant advantages compared to conventional finite element 
approximations. In particular, a widely cited strength of such methods is 
improved accuracy for the same computational cost. Yet, few attempts 
have been made to directly assess the quality of the results obtained dur- 
^.^ ing the simulation by evaluating an estimate of the discretization error. 

,^^ Here we propose a recovery type error estimator based on an enhanced 

y^ recovery technique. The salient features of the recovery are: enforcement 

^H of local equilibrium and, for singular problems a "smooth+singular" de- 

C/3 composition of the recovered stress. We evaluate the proposed estimator 

O on a number of test cases from linear elastic structural mechanics and ob- 

tain precise error estimations whose effectivities, both at local and global 
,_^ levels, are improved compared to recovery procedures not implementing 

^ these features. 

oo 

1^^ Keywords: smoothed finite element method; error estimation; statical admissibility; SPR- 

^S| CX; singularity; recovery 
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The smoothed finite element for mechanics problems was introduced in 2006 
by Liu et al., [1]. The main idea of the method is to relax the compatibility 
condition at the element level by replacing the standard compatible strain by 
its smoothed counterpart. The smoothing operation can be performed over do- 



^ mains of various shapes which can be obtained by dividing the computational 

domain into non- overlapping smoothing domains. These domains can be ob- 
tained by subdividing the elements (cell-based smoothing) as in [1, 2, 3, 4, 5], or 
using edge [6, 7] or node-based geometrical information [8]. Each method has 
several advantages and drawbacks, summarized in, e.g. [9, 10], but the strongest 
motivation for smoothed finite elements is certainly revealed in its enhancement 
of low order simplex elements (e.g. linear triangles and tetrahedral), alleviating 
overstifFness, locking and improving their accuracy significantly [11, 12, 6]. 
The applications of strain smoothing in finite elements are wide. Since the 
introduction of the smoothed finite element method (SEEM), the convergence. 
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the stability, the accuracy and the computational complexity of this method 
was studied in [2, 3] and the method was extended to treat various problems 
in solid mechanics such as plates [13], shells [14] and nearly incompressible 
elasticity [12]. Recently, Bordas et al. [15, 16] combined strain smoothing with 
the XFEM to obtain the Smoothed eXtended Finite Element Method to solve 
problems with strong and weak discontinuities in 2D continuum. 
In this paper, we focus on the cell-based smoothing, a review of which is pro- 
vided in [4, 12, 17] along with applications to plates, shells, three dimensional 
continuum and a coupling with the extended finite element method with appli- 
cations to linear elastic fracture (continuum, plate). 

The development of new numerical techniques based on the finite element method 
(e.g. GFEM, XFEM, ...) aims at obtaining more accurate solutions for engineer- 
ing problems. Despite the accuracy gained with the new techniques, numerical 
errors, especially the discretization error, are always present and have to be 
evaluated. Consolidated accuracy assessment techniques developed in the FE 
framework are commonly adapted to the framework of these new techniques 
[18, 19, 20, 21, 22]. As in any numerical method, the smoothed FEM approxi- 
mation introduces an error that needs to be controlled to guarantee the quality 
of the numerical simulations. Although an adaptive node-based smoothed FEM 
has been developed in [23], a rather simple error estimator using a recovery 
procedure which initially is only valid for NS-FEM is used to guide the adaptive 
process. The technique evaluates a first-order recovered strain field interpolating 
the nodal values by means of the linear FEM shape functions. 
The urge for quality assessment tools for smoothed FEM approximations, the 
promising results in [23] obtained with a rather simple technique and the ex- 
perience of the authors in the development of high quality recovery-based error 
estimators in the FEM and XFEM contexts, motivates this paper where we esti- 
mate, a posteriori, the approximation error committed by cell-based smoothed 
finite elements. The method used for a posteriori error estimation relies on 
the Zienkiewicz and Zhu error estimator [24] commonly used in FEM, together 
with a recovery technique recently developed by the authors, specially tailored 
to the analysis of enriched approximations containing a smooth and a singular 
part and which locally enforces the fulfilment of equilibrium equations. The 
technique known as SPR-CX [21, 25] was shown to lead to very good effectivity 
indices in FEM and XFEM. 

The paper is organised as follows: In Section 2, the boundary value problem 
of linear elasticity is briefly introduced and the approximate solution using the 
smoothed finite element method is presented. In Section 3, we discuss basic con- 
cepts related to error estimation, in especial recovery-based techniques. Section 
4 is devoted to the proposed enhanced recovery technique and its application to 
SFEM approximations. Numerical examples are presented in Section 5 and the 
main concluding remarks in Section 6. 

2 Problem statement and SFEM solution 

2.1 Problem statement 

Let us consider the 2D linear elasticity problem. The unknown displacement 
field u, taking values in 51 C M^, is the solution of the boundary value problem 



given by 

-V-a(u)=b infl (1) 

o- (u) • n = t on Tn (2) 

u = on Yd (3) 

where r^r and F^j denote the Neumann and Dirichlet boundaries with dfl = 
Tn U Td and Tn H Tjj — 0. The Dirichlet boundary condition in (3) is assumed 
to be homogeneous for the sake of simphcity. 
The weak form of the problem reads: Find u G V such that 

VveV^ a(u,v) = Z(v), (4) 

where V is the standard test space for the elasticity problem such that V — 
{v|vG[//i(r!)]2,v|r„(x) = 0}, and 

a(u,v) := / eiuf'De{w)dn = f (t(u)'^D- V(v)df} (5) 

Jn Jn 

l{-v):= I h^wdVt+ [ t^vdF, (6) 



where D is the elasticity matrix of the constitutive relation cr = De , cr and s 
denote the stress and strain operators. 

2.1.1 Singular problem 

Figure 1 shows a portion of an elastic body with a reentrant corner (or V-notch) , 
subjected to tractions on remote boundaries. No body loads are applied. For 
this kind of problem, the stress field exhibits a singular behaviour at the notch 
vertex. 



Figure 1: Sharp reentrant corner in an infinite half-space. 

The analytical solution corresponding to the stress distribution in the vicinity 
of the singular point is a linear combination of the singular and the non-singular 



terms. It is often claimed that the term with a highest order of singularity dom- 
inates over the other terms in a zone surrounding the singular point sufficiently 
closely. The analytical solution to this singular elastic problem in the vicinity 
of the singular point was first given by [26] and is described, for example, in 
[27, 28]. Here, we reproduce those expressions for completeness. In accordance 
with the polar coordinate system of Figure 1, the displacement and the stress 
fields at points sufficiently close to the corner can be described as: 

u(r, <j>) = i^ir^^*i(Ai, (/)) + i^„r^"*„(A„, 0) (7) 

cr(r, 0) = i^iAir^i-i$i(Ai, cb) + i^iiAiir^"-i*ii(A„, 0) (8) 

where r is the radial distance to the corner. Am (with m = I, II) are the eigen- 
values that determine the order of the singularity, ^m and $„ are a set of 
trigonometric functions that depend on the angular position (/), and Km are 
the so-called generalised stress intensity factors (GSIFs). The generalized stress 
intensity factor (GSIF) is a multiplicative constant that depends on the loading 
of the problem and linearly determines the intensity of the displacement and 
the stress fields in the vicinity of the singular point. Hence, the eigenvalues A 
and the GSIFs K define the singular field. 

The eigenvalue A is easily known because it depends only on the corner angle a, 
and can be obtained as the smallest positive root of the following characteristic 
equations [26]: 

sin AjQ! + Ai sina = (9) 

sin AiiQ! — All sina = (10) 

The smallest positive root yields the highest order of singularity and determines 
the term that dominates the elastic fields given by (7) in the vicinity of the notch 
vertex. Strictly speaking, (9) corresponds to the symmetric part of the elastic 
fields with respect to = (i.e. the bisector line BB in Figure 1) and (10) to 
the antisymmetric solution. These solutions are also called mode I and mode II 
solutions, respectively. The trigonometric functions for the mode I displacement 
and stress fields in (7, 8) are given by [27]: 



*i(Ai 



*I,a;(Al,0)' 
*I,y(Al,</'), 

= J_ /(k - g(Ai + 1)) COS Ai(/)-Aicos(Ai- 2)0)1 , . 

2m I (« + 0(Ai + 1)) sin Ai<^ + Ai sin(Ai - 2)0 J ^^^' 

{*La;a;(Al,0)' 
*i,to(Ai,0) 

' (2 - Q(Ai + 1)) cos(Ai - 1)0 - (Ai - 1) cos(Ai - 3)0] 
(2 + Q(Ai + 1)) cos(Ai - 1)0 + (Ai - 1) cos(Ai - 3)0 \ (12) 
g(Ai + 1) sin(Ai - 1)0 + (Ai - 1) sin(Ai - 3)0 J 

where n is the shear modulus and k is the Kolosov constant, defined as functions 
of E (Young's modulus) and v (Poisson's coefficient) according to the following 
expressions: 

J7I [3 — 4;/ plane strain 

M — W7 ^ ' K ^ i 3 — V 

2[1 + ly) plane stress 

Kl + i^ 



In the same way, for mode II we have: 

^ J_ f (K-(5(Aii + l))sinAii0- Aiisin(An-2)(/) 1 , , 

2m !-('« + 0(Aii + l)) cos A„0- All cos(Aii-2)(/)/ ^ ' 

{*II,a;a;(AlI,(/))' 
*II,a;j,(AlI,0)^ 

' (2 - 0(Aii + 1)) sin(Aii - 1)0 - (An - 1) sin(Aii - 3)(j)] 
(2 + 0(Aii + 1)) sin(Aii - 1)0 + (An - 1) sin(Aii - 3)0 \ (14) 
-(3(Aii + 1) cos(Aii - 1)0 + (An - 1) cos(Aii - 3)0 J 

Note that the components of the displacement and the stress fields are expressed 
in Cartesian coordinates. In addition, Q is a constant for a given notch angle: 

cos((Ai-l)|) sin((A„-l)|) 

Qi 7 ^Y ' '3" 7 a\ (1^) 

cos((Ai + l)-j sin(^(An + l)-J 

2.2 FEM solution with strain smoothing 

2.2.1 Finite element formulation 

Let u'' be a finite element approximation to u. The solution lies in a func- 
tional space V'^ C V associated with a mesh of isoparametric finite elements of 
characteristic size h, and it is such that 

Vv'' e y'' a(u'',v'') = ?(v'') (16) 

Using a variational formulation of the problem in (1-3) and a finite clement 
approximation u'' — Nu"^, where N denotes the shape functions of order p, we 
obtain a system of linear equations to solve the displacements at nodes u*^ : 

KU = f (17) 

where K is the stiffness matrix, U is the vector of nodal displacements and f is 
the load vector. 

2.2.2 Strain smoothing in FEM 

Inspired by the work of Chen et at, in 2001 [29] on stabilized conforming 
nodal integration (SCNI), Liu et at, [1] introduced the smoothed finite ele- 
ment method. The idea behind SCNI/SFEM is to write a strain measure as a 
spatial average of the standard strain field. To do so, the elements are divided 
into smoothing cells over which the strain is smoothed. By the divergence theo- 
rem, integration over the element is transformed to contour integration around 
the boundaries of the subcell. A particular element can have a certain num- 
ber of smoothing cells and depending on that number, the formulation offers 
a range of different properties [1, 2, 3]. Interested readers are referred to the 
literature [2, 3, 12] for a detailed description of the method, its variants [7, 17, 8] 



and its convergence properties [3, 11]. The strain field e^^-, used to evaluate the 
stiffness matrix is computed by a weighted average of the standard strain field 
e^-. At a point xc in an element fi'*, 



S% (xc ) = / e% (x) $ (x - xc ) dx 



(18) 



where $ is a smoothing function that generally satisfies the following proper- 
ties [301 



$ > and / $(x)rfx = 1 

One possible choice of $ is given by: 
1 



$ = — — in Vic and $ = elsewhere 
Ac 



(19) 



(20) 



where Ac is the area of the subcell. To use (18), the subcell containing point 
X(7 must first be located in order to compute the correct value of the weight 
function $. 

The discretized strain field is obtained through the smoothed discretized gradi- 
ent operator or the smoothed strain displacement operator, Be, defined by 



£'*(xc.)-Bc(xc)q 



(21) 



where q are the unknown displacements coefficients defined at the nodes of the 
finite element, as usual. The smoothed element stiffness matrix for element e is 
computed by the sum of the contributions of the subcells^ 



K*^ 



nc 



C=l 



Bj^DBcdQ 



/ I ^c^"ca" = 2^ B^DBc 



c=i 



« nc 

/ dn^Y^^c^^cAc (22) 

J ^c C—1 



where nc is the number of the smoothing cells of the element. The strain 
displacement matrix Be is constant over each flc and is of the following form 



Br 



Bci B 



C2 



B 



C3 



B 



C4 



(23) 



where for all shape functions / e {1, . . • ,4}, the 3x2 submatrix Bp/ repre- 
sents the contribution to the strain displacement matrix associated with shape 
function / and cell C and writes 



VJ G {1, 2, . . . , 4}, VC G {1, 2, . . . nc}Bci = 



Sc 



nx 
n. 



y 



{x)Ni{x)dS 



(24) 

or, since (24) is computed on the boundary of ilc and one Gaufi point is sufficient 
for an exact integration (in the case of a bilinear approximation): 



tThc subccUs Qc form a partition of the element Q^. 



/ 



^ nb 



Nj (xf ) n, 

Nj (xG) n, 

TV, (xG) n^ TV, (xG) n, ^ 



?,^ (25) 



where nj, is number of edges of the subceU, {nx,ny) is the outward normal to 

cf and l^ 
length of r^, respectively. 



the smoothing cell, i^c, ^f and l^ are the center point (Gaufi point) and the 



3 Error estimation in the energy norm 

The discretization error in the standard finite clement approximation is defined 
as the difference between the exact solution u and the finite element solution 
u'* : e = u — u'' . Since the exact solution is in practice unknown, in general, 
the exact error can only be estimated. To obtain an estimation of e, norms 
that allow a better global interpretation of the error arc normally used. The 
Zienkiewicz-Zhu [24] error estimator defined as 

{(T* - a'^f B-^ {(T* - (t'') dn (26) 

n 

relies on the recovery of an improved stress field <t*, which is supposed to be 
more accurate than the FE solution cr^, to obtain an estimation of the error in 
energy norm ||ees||. The domain fl could refer to the full domain of the problem 
or a local subdomain (element). 

The recovered stress field cr* is usually interpolated in each element using the 
shape functions N of the underlying FE approximation and the values of the 
recovered stress field calculated at the nodes cr* , given by: 

<T*(x)=£iV,(x),TKxz), (27) 

where rig is the number of nodes in the element under consideration and <Tj(x/) 
are the stresses provided by a recovery technique at node /. The supercon- 
vergent patch recovery technique (SPR) proposed by Zienkiewicz and Zhu [3f] 
is commonly used to evaluate the components (j — xx, yy, xy) of cr| using a 
polynomial expansion, cr/j = pa. This expansion is defined over a set of con- 
tiguous elements connected to node / called patch, where p is the polynomial 
basis and a are the unknown coefficients obtained using a least squares fitting to 
the values of the FE stresses evaluated at integration points in the patch, being 
p, normally, of the same order as the interpolation of displacements. The ZZ 
error estimator is asymptotically exact (i.e. the approximate error converges to 
the exact error as the mesh size goes to zero) if the recovered solution used in 
the error estimation converges at a higher rate than the finite element solution 
[31, 32]. As it can be seen in (26), the accuracy of the error estimate is closely 
related to the quality of the recovered field. For this reason, several techniques 
have been developed aiming to improve the quality of cr* . The authors have 



proposed different techniques mostly for the FEM/XFEM context as the ex- 
tended moving least squares recovery (XMLS) and the extended global recovery 
techniques proposed by Duflot and Bordas in a series of papers [19, 33, 34], the 
SPR-C and the SPR-CX by Rodenas et al. [35, 21], which were used later as the 
basis for the development of recovery-based error bounding techniques [25, 36]. 
The next section presents the SPR-CX technique which improves the recovered 
field by enforcing equilibrium and effectively dealing with singular fields. 
Remark: In mathematics is common to consider that one can only speak about 
an error estimator if sharp or at least approximated upper - and desired - 
also lower error bounds can be proven, reserving the word indicator when the 
technique does not necessarily bound the error. However, this terminology is 
not general and many other authors, usually from the engineering community, 
use the term error estimator even when the technique is not able to provide 
error bounds. This is the case for example in [24, 37, 38] and also our case. 

4 SPR-CX recovery technique 

The SPR-CX recovery technique first introduced by Rodenas et al. [21] is an 
enhancement of the superconvergent patch recovery technique in [31], which 
incorporates the ideas proposed in [35] to guarantee the exact satisfaction of 
the equilibrium locally on patches. In [21, 25] a set of key ideas are proposed 
to modify the standard SPR allowing its use with singular problems. 
The recovered stresses cr* are directly evaluated at a sampling point (e.g an 
integration point) x through the use of a partition of unity procedure, prop- 
erly weighting the stress interpolation polynomials obtained from the different 
patches formed at the vertex nodes of the element containing x: 

^*(x)=£7V,(x)^Kx), (28) 

7=1 

where Ni are the shape functions associated to the vertex nodes n^. To obtain 
the nodal values cri^ we solve a least squares approximation of the stresses eval- 
uated at a set of sampling points distributed within the domain of the patch 
of node / (elements connected to /). In FEM, such points usually correspond 
to the integrations points used in the finite clement approximation. In SFEM, 
we map the constant strains at each subcell to a 2 x 2 Gauss quadrature dis- 
tribution in the subcell used as sampling points. This way we have a sufficient 
number of points at each patch to solve the linear system of the least squares 
approximation, see Figure 2. Note that as in the other versions of the SFEM 
(NS-FEM, ES-FEM) the elements are also subdivided into subcells, a similar 
approach can be used to perform the mapping of the stresses to sampling points. 
Therefore, the proposed error estimation technique can be used with all SFEM 
implementations. 

One major modification of the original SPR technique is the introduction of a 
splitting procedure to perform the recovery. For singular problems the exact 
stress field cr is decomposed into two stress fields, a smooth field cTsiyio and a 
singular field cTsing'- 

^ ^ s7no I ^ sing- V ^j 



X X 11 X X 

* !' * 

X X i! X X 



j Subcell 
• Subcell centre 
X Sampling points 



Figure 2: Distribution of stress sampling points at each subcell in a 2-subcell 
quadrilateral element used in the stress projection. 



Then, the recovered field cr* required to compute the error estimate given in 
(26) can be expressed as the contribution of two recovered stress fields, that is, 
smooth (r\.„^a ^i^*^ singular cr'^^^g (see Figure 3): 



(30) 
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Figure 3: Evaluation of the recovered field at different patches. 

For the recovery of the singular part, the expressions in (8) which describe the 
asymptotic fields in the vicinity of the singular point are used. The radius p that 



defines the splitting has to be large enough as to ensure that in the first mesh 
at least the patch of elements corresponding to the node located at the singular 
point is within the singular + smooth decomposition area. To evaluate cr*ginq 
from (8) wc first obtain estimated values of the stress intensity factors K\ and 
Kii using the interaction integral as indicated in [21, 25, 39]. The recovered part 
cr*j„ is an equilibrated field as it satisfies the internal equilibrium equations. 
Once the field cr*^^ has been evaluated, an FE approximation to the smooth 
part cr^smo can be obtained subtracting cr^^^ from the raw FE field: 

^srno = "" ^ ^ sing- \^^) 

Then, the field c*^^ is evaluated by applying an SPR-C recovery procedure 
over the field cr^rno- 

In order to obtain an equilibrated recovered stress field, the SPR-CX enforces 
the fulfilment of the equilibrium equations locally on each patch. The constraint 
equations arc introduced via Lagrange multipliers into the linear system used to 
solve for the coefficients of the polynomial expansion of the recovered stresses 
on each patch. These include the satisfaction of the: 

• Internal equilibrium equations. 

• Boundary equilibrium equation: A point collocation approach is used to 
impose the satisfaction of a second order approximation to the tractions 
along the Neumann boundary. 

• Compatibility equation: This additional constraint is also imposed to fur- 
ther increase the accuracy of the recovered stress field. 

To evaluate the recovered field, quadratic polynomials are used in the patches 
along the boundary and crack faces, and linear polynomials for the remaining 
patches. As more information about the solution is available along the boundary, 
polynomials one degree higher are useful to improve the quality of the recovered 
stress field. 

The enforcement of equilibrium equations provides an equilibrated recovered 
stress field locally on patches. However, the process used to obtain a continuous 
field a* shown in (28) introduces a small lack of equilibrium as explained in 
[25, 36]. The reader is referred to [25, 36, 40, 41] for details. 
In order to estimate the error in SEEM approximations we can follow a similar 
procedure. To build the patches we use the topological information of the SEEM 
discretization. The recovered stress field is evaluated at the centre of the subcells 
and then projected to the sampling points as explained before. 
Remark: The recovery method proposed in this paper is general, and could 
also be applied, although this is beyond the scope of this paper, to problems with 
corner singularities at triple junctions in polycrystalline materials made up of 
orthotropic grains to estimate the error of extended finite element formulations 
such as those recently proposed in [42, 43]. 

5 Numerical results 

In this section, numerical tests considering 2D benchmark problems with exact 
solution have been used to investigate the quality of the proposed technique. 
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The performance of the technique has been evaluated using the effectivity index 
of the error in energy norm, both at global and local levels. Globally, we have 
considered the value of the effectivity index 6 given by 



(32) 



where ||e|| denotes the exact error in the energy norm, and ||eg^|| represents 
the evaluated error estimate. At element level, the distribution of the local 
effectivity index D, its mean value rndDl) and standard deviation a{D) have 
been analysed, as described in [35]: 



= 9" -1 if e^' >1 
D = l- — if 0" <1 



with 



1^ 
lle'^ 



(33) 



Note that 6''^ e (0, 1) when the error is underestimated and 9'^ G (1, +oo) when 
it is overestimated. The definition of D fairly compares the underestimation of 
the error {D < 0) and the overestimation (D > 0). The good local behaviour of 
the estimates results in values of D close to zero. The global effectivity index 
is used to evaluate global results. The mean value to(|-D|) and the standard 
deviation (t{D) of the local effectivity are also used to evaluate the global quality 
of the error estimator as these parameters are useful to take into account error 
compensations in the evaluation of 0. 

5.1 Thick-wall cylinder subjected to an internal pressure. 

The geometrical model for this problem is shown in Figure 4. Due to symmetry 
conditions, only one part of the section is modelled. Plane strain conditions are 
assumed. 




Plane strain 

a=l 

b=5 

E=3.10' 
v=0.3 



Figure 4: Thick-wall cylinder subjected to internal pressure. 



The exact solution to this problem is given by the following expressions, where 
for a point {x,y), c — b/a, r = yjx^ + y"^ and </> = arctan(y/x) the radial 
displacement is given by 



i;^(^<--'4 



(34) 
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Stresses in cylindrical and cartesian coordinates are 



P 

C2-1 
P 



CTt = 



-1 



1- 



1 



62 

?! 



CTxx = cTr cos2((/)) + cTt sin^(0) 
(jyj, = (Jr sin^(0) + <7t cos^{(j)) 

(Txy = (o-r - (Tt) Sm{(j)) COs{(j)) 



(35) 



c^ - 1 

A sequence of uniformly refined meshes of linear quadrilateral elements have 
been used for the analyses. The material parameters are Young's modulus 
E = 'ix 10^ and Poisson's ratio v = 0.3. In the case of the SFEM approximation, 
the element is divided into 4 subcells. However, the influence of the number of 
subcells on the global/local error level is also studied in a later analysis. Figure 
5 shows the exact error for the raw SFEM and the recovered stress fields for the 
three stress components and the von Mises stress. It can be seen that the error 
in the recovered field is significantly smaller than the error for the raw stress 
solution. 
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-1000 1000 



Figure 5: Cylinder under internal pressure. Exact stress error for the SFEM 
and recovered fields considering the three stress components and the von Mises 
stress. 

Figure 6 shows the distribution of the local effectivity index for a sequence of 
uniformly refined meshes. The local effectivity values are within a very narrow 
range and they improve with mesh refinement. The distribution of the D is 
homogeneous and good results are obtained along the boundary. 
Figure 7 shows the convergence of the estimated error in energy norm using 
two different configurations of the recovery procedure: the curve SPR-CX for 
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Figure 6: Cylinder under internal pressure. Distribution of the effectivity index 
D. 



an equilibrated recovered field and the curve SPR for a non-equilibrated recov- 
ery, the exact error (with a convergence rate s = 0.5, which equals the optimal 
convergence rate in the energy norm with respect to the number of degrees of 
freedom for a smooth problem) is shown for comparison. We have solved the 
same problem using a standard FEM approximation and an equilibrated recov- 
ery technique (curve SPR-CX (FEM)) to estimate the error in that solution, 
the results are also included in the figure. The FEM values using the SPR-CX 
recovery converge with a rate of 0.49, while for the SFEM using SPR-CX and 
standard SPR show an average convergence rate of 0.49 and 0.42, respectively. 
The equilibrated SFEM and FEM error estimates obtained using the SPR-CX 
technique are both very close to their corresponding exact errors and converge 
with the same rate. If we do not consider equilibrium constraints during the 
recovery (SPR curve) the error is underestimated and the convergence rate is 
lower. These results clearly show the importance of the use of equilibrated 
recovery techniques for accurate error estimations. 




- SPR-CX 
-SPR 

- SPR-CX (FEM) 
-exact, (FEM) 
-exact, (SFEM) 



Figure 7: Cylinder under internal pressure. Convergence of the estimated error 
llee^ll for the SFEM using 4 subcells (SPR-CX, s = 0.49), the SFEM without 
enforcing equilibrium (SPR, s = 0.42) and the FEM (SPR-CX (FEM), s = 
0.49). The exact error is shown for comparison. 
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Figure 8: Cylinder under internal pressure. Evolution of the estimated error 
||ees|| and convergence rate s for the SFEM using 4 subcells (SPR-CX, Savg — 



0.49) and for the error in the recovered solution (||e*||, 



^avg 



0.65). 



In Figure 8 we represent the evolution of the estimated error for the SFEM 
solution using the enhanced recovery (curve SPR-CX) and the exact error of 
the recovered field || e* || = || u— u* || . The error in the recovered field has a higher 
convergence rate, which is in agreement with the expected quality for the field 
(T* . According to [32] , this also serves to verify the asymptotic exactness of the 
proposed error estimator. 

Figure 9 shows the evolution of global indicators 9, m{\D\) and (y{D) for the 
SFEM using equilibrated recovery (SPR-CX curve), without equilibrium con- 
straints (SPR curve) and the standard FEM with equilibrium (SPR-CX (FEM) 
curve). The equilibrated SFEM and the FEM recoveries exhibit similar results, 
with good effectivity of the error estimator and decrease of to(|Z?|) and (^{D) for 
finer meshes. The non-equilibrated SFEM recovery (SPR curve) shows worse 
results and converges with a slower rate. 

The use of different numbers of subcells for the SFEM approximation is also 
considered for comparison. Figure 10 shows the convergence of the estimated 
error in energy norm for two, four and eight subcells. All the curves exhibit the 
same convergence rate (s = 0.49), close to the theoretical value s = 0.5 
Figure 11 shows the evolution of global indicators 9, m{\D\) and cr(D) for two, 
four and eight subcells. The effectivity indices for all the subcells types shown 
converge asymptotically to the theoretical value and are very accurate (1.08 > 
9 > 1). The local effectivity index goes to zero at the same rate as shown in the 
curves 7n(|I?|) and cr{D). 

In Figure 12 we show the influence of the order of the polynomial expansion 
used for the local recovery on patches. We compare the evolution of the global 
parameters for first order polynomials, previously represented in Figure 11, with 
the corresponding curves considering second order polynomials. We can see that 
the increase of the polynomial order does not produce and improvement of the 
effectivity. Local behaviour in 7ti(|I?|) and (y{D) indicate even worse results as 



14 



10- 



10- 



1.1 



0.9 



0.8 



0.7 




m(|D|) 



10^ 





SPR-CX 

SPR 

SPR-CX (FEM) 



10=* 



dof 



aiD) 





10- 



- SPR-CX 
-SPR 

- SPR-CX (FEM) 



dof 



Figure 9: Cylinder under internal pressure. Global indicators 9, m{\D\) and 
(t{D) for SPR-CX, SPR and SPR-CX (FEM). 




Figure 10: Cylinder under internal pressure. Convergence of the estimated error 
||ees|| for elements with two (2SC), four (4SC) and eight (8SC) subcells. 



we increase the number of degrees of freedom. This is in correspondence with 
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Figure 11: Cylinder under internal pressure. Global indicators 0, to(|_D|) and 
a{D) for elements with two (2SC), four (4SC) and eight (8SC) subcells. 



previous results observed in the FEM context [44], where an increase of the 
polynomial order not necessarily derived in better effectivities. 

5.2 L-shape domain under mode I load. 

This singular problem consists of a portion of an infinite L-shaped domain. The 
model is loaded on the boundary with the tractions corresponding to the first 
symmetric term of the asymptotic expansion that describes the exact solution 
under mode I loading conditions around the singular vertex, see Figure 13. The 
exact values of boundary tractions on the emphasized boundaries have been 
imposed in the FE analyses. 

The exact displacement and stress fields for this problem are given by (7, 8). For 
a = 37r/2 one obtains Ai = 0.544483736782464, An = 0.908529189846099 and 
Q = 0.543075578836737. Exact values of the GSIFs have been taken as Xj = 1 
and A'n = 0. The material parameters are Young's modulus E — 1000, and 
Poisson's ratio v = 0.3. The splitting radius is p = 0.5. To evaluate the 
SIF we use the equivalent interaction integral defined with a plateau function 
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Figure 12: Cylinder under internal pressure. Global indicators 6*, to(|D|) and 
(j{D) for elements with two (2SC), four (4SC) and eight (8SC) subcells using 
1st and 2nd order polynomials for the recovery on patches. 



with radius 0.9 as indicated in [21]. As the analytical solution of this problem 
is singular at the reentrant corner of the plate, we apply the singular-fsmooth 
decomposition of the stress field as explained in Section 4. We use the expression 
in (8) and the estimated GSIF evaluated using the interaction integral to obtain 
an approximation of the singular part Using ■ 

Figure 14 shows the distribution of the local effectivity index for the sequence of 
graded meshes in Figure 15. In the figure, the local index D decreases with the 
refinement of the meshes and the obtained values are within a narrow range. The 
first mesh of the sequence is an uniform mesh. This kind of meshes is known to 
produce pollution error for singular problems. The highly underestimated areas 
at the right of the plate are explained by this pollution error, which cannot be 
controlled with local techniques. 

Figure 16 shows the convergence of the estimated error in energy norm for 
different configurations of the recovery procedure: SPR-CX that consider equi- 
librium and stress decomposition, SPR-X that considers stress decomposition 
only, SPR-C that consider equilibrium only, and a conventional SPR. The exact 
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Figure 13: L-shaped domain under mode I load. 




-0.2 0.2 



Figure 14: L-shaped domain under mode I load. Distribution of the efFectivity 
index D. 




Figure 15: L-shaped domain under mode I load. Sequence of graded meshes. 



error is shown for comparison. 

Figure 17 shows the evolution of global indicators 9, m{\D\) and cr(-D) for the 
different configurations: SPR-CX, SPR-X, SPR-C, SPR. The best resuhs are 
for the SPR-CX. The SPR-C and the SPR cannot accurately recover the field 
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Figure 16: L-shaped domain under mode I load. Convergence of the estimated 
error ||ees|| for for different configurations of the recovery technique: SPR-CX, 
SPR-X, SPR-C, SPR. The exact error is shown for comparison. 



close to the singularity. The SPR seems to provide good global effectivity re- 
sults, however, this is only due to compensation between underestimated and 
overestimated areas. The real behaviour is clear when analysing the evolution 
of ?7i(|D|) and (j{D) and even more patent if we represent D as seen in Figure 
18, where we represent the results for the different configurations of the recovery 
technique. The effect of equilibrium enforcement and singular decomposition is 
clearly shown. A more homogeneous distribution of D is obtained when using 
SPR-CX. SPR and SPR-X are not equilibrated along the boundary and there- 
fore they underestimate the error in the elements along it. SPR and SPR-C 
overestimates the error in the vicinity of the reentrant corner. 

6 Conclusions 

In this paper, an a posteriori recovery-based error estimator which makes use 
of a modified version of the SPR technique previously used in the FEM and 
XFEM contexts has been adapted to the smoothed finite element method. The 
recovery technique considers the local enforcement of equilibrium equations and 
the singular + smooth decomposition of the stress field for singular problems. 
The technique has been applied to the cell-based smoothed FEM but could also 
be used with the node-based and edge-based smoothed FEM implementations. 
The numerical results presented in this paper show that the method yields ac- 
curate estimations of the error in the energy norm both locally and globally. 
It can be inferred that enforcing equilibrium constraint is required to obtain 
accurate results. The influence in the recovery of the number of subcells used to 
formulate the smoothed finite elements is also studied, showing that the tech- 
nique performs adequately for the different configurations. The error estimator 
based on the use of the SPR-CX recovery technique accurately captures the 
discretization error both in smooth and singular problems. Moreover, it could 
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Figure 17: L-shaped domain under mode I load. Global indicators 6', rndDj) and 
cr{D) for different configurations of the recovery technique: SPR-CX, SPR-X, 
SPR-C, SPR. 



be used to guide /i-adaptive refinements and the recovered field <t* can be used 
as an enhanced solution, more accurate than the stress field provided by the 
approximation. 

Future work includes the comparison of the proposed estimators with other re- 
cently developed error estimators for the extended finite element method when 
dealing with singular problems [19, 33, 34] for three dimensional fracture prob- 
lems, the focus of our current research. We will also analyse the behaviour of 
strain smoothing for real-life three dimensional fracture mechanics problems, 
which is the topic of the EPSRC project which funded this work. 
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